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1.  INTRODUCTION 

Consider  the  model 
t0 

y(n)  =  l  a.  exp(io).n)  +  w(n),  n«l,2,  ...,N  (1.1) 

k=l  K  K 

where  f  0  are  unknown  amplitudes,  e  (0,2it),  (k  =  1,2,  ...,t»),  are 

unknown  frequencies  and  i  =  /-].  Usually  we  assume  that  the  noise  w(n)  has 

2 

mean  zero  and  finite  variance  a  .  The  above  model  is  of  interest  in  tne 
area  of  signal  processing.  Under  the  above  model,  it  is  of  interest  to 
estimate  the  unknown  parameters.  Even  when  tp  is  known,  it  is  difficult 
to  find  the  least  square  estimates  of  a^'s  and  u^'s  since  it  would  involve 
solving  a  system  of  nonlinear  equations  with  exponential  functions.  To 
avoid  this  difficulty,  various  methods  have  been  developed  in  the  literature, 
such  as  linear  prediction  (LP)  method,  those  based  upon  using  principal 
eigenvectors  of  estimated  cross-spectral  density  matrices  (Liggett  (1973)) 
and  the  forward-backward  linear  prediction  (FBLP)  method  (Nuttall  (1976) 
and  Ulrych  and  Clayton  (1976)). 

In  the  LP  procedure,  we  have  still  to  solve  a  polynomial  equation  whose 
degree  would  be  rather  high  although  we  do  not  solve  a  system  of  exponential 
equations.  Also,  it  does  not  work  well  for  the  case  of  low  signal-to-noise 
ratio  (SNR).  Tufts  and  Kumaresan  (1982)  made  modification  to  the  original 
FBLP  method,  and  showed  by  simulation  that  the  modified  FBLP  works 
much  better  than  the  original  FBLP  when  the  SNR  is  relatively  low.  However, 
it  still  involves  solving  a  high  degree  polynomial  equation. 

In  the  present  paper,  we  investigate  the  estimation  of  both  the 
number  of  signals  and  the  amplitudes  and  frequencies  of  the  signals,  and 
study  a  method  which  we  refer  to  as  equivariation  linear  prediction 


S 


■  V.V.VJ 


(EVLP).  Using  this  approach,  we  can  at  the  same  time  find  the  estimates  of 
number  of  signals  and  the  frequencies  and  then  using  the  usual  least  square 
method  to  get  the  estimates  of  the  amplitudes.  In  this  method,  we  need  only 
to  solve  a  polynomial  equation  with  the  lowest  degree.  In  Section  2  we 
shall  state  this  method.  In  Section  3  we  shall  prove  the  strong  consistency 
of  this  procedure.  In  Section  4  we  give  the  limiting  distribution  of  vari¬ 
ous  estimates,  given  in  previous  sections.  In  Section  5  we  qive  a  further 
discussion  on  the  FBLP  and  the  modified  one.  In  this  section  we  will  show 
the  reason  why  (theoretically)  the  modified  FBLP  works  better  than  the 
original  one  when  the  SNR  is  relatively  low  and  the  sample  size  is  small, 
and  show  that  these  two  methods  will  become  equivalent  when  the  sample  size 
goes  to  infinity.  More  importantly,  we  shall  point  out  that  both  these  two 
methods  do  not  provide  consistant  estimation  of  the  frequencies,  and  we 
shall  pose  a  further  modification  on  FBLP  such  that  the  procedure  is  con- 
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2.  DETERMINATION  OF  THE  NUMBER  OF  THE  FREQUENCIES 
AND  ESTIMATION  OF  THE  FREQUENCY  PARAMETERS 

Suppose  that  the  data  sequence  y(n)  is  given  by  the  formulas 
t0 

y(n)  =  l  a,  exp(ioj.n)  +  w(n) ,  n  =  l,2,...,N,  (2.1) 

j=l  J  J 

where  i  =  /T,  {a.}  is  a  set  of  unknown  complex  amplitudes,  {w- }  is  a  set 
J  J 

of  unknown  angular  frequencies,  and  (w(n)}  is  a  sequence  of  i.i.d.  complex 
random  noise  variables  such  that 

Ew ( 1 )  =  0 ,  Ew(1)wTTT  =  o2,  E | w(  1 )  | 4  <  «>  (2.2) 

p 

with  o  unknown.  We  assume  that  w.  f  w.  if  j  f  k,  and  u>.  e  (0,2tt)  for  all  j. 

J  K  J 

In  this  paper,  we  are  primarily  interested  in  determining  tQ  and  esti¬ 
mating  the  frequency  parameters  c u..  Once  w.'s  are  accurately  determined, 

J  J 

the  a-'s  can  be  found  by  a  linear  least  squares  fit  to  the  data. 

w 

To  determine  tg,  the  true  number  of  different  angular  frequencies,  it 
is  prescribed  a  prior  that  tg  _<  T  <  °®. 

For  any  nonnegative  integer  t  <  T,  write  complex  vector  b^  as 


b(t)  =  (bj,t},  ...,b[t))‘. 


(2.3) 


Put1 


St  =min{^  l  |  l  b^t)y(n-0!2:  ||b(t)||=l}. 


n=t+l  l=0 


t  =  0,1,2 . T, 


(2.4) 


★ 


Rao  (1986)  also  remarked  about  finding  S^.,  for  given  t,  to  estimate  the 
frequencies. 


where  ||h  t^||  =  (  I  |b^  |2)ly/2.  Take  CN  satisfying  the  following  condi¬ 


tions  : 


Tim  CN  =  0  and  lim/N  CN//log  log  N  = 
N-*=°  N-**>  '  ~ 


(2.5) 


Then  we  can  find  a  nonnegative  integer  t  s  ^  <  T  which  minimize 


Rt  "  St  +  tCN*  1  ‘  °»  1  ’  •••»  T» 


(2.6) 


and  use  tQ  as  an  estimate  of  tQ. 


Further,  we  can  find  a  unit  (tg+l)xl  complex  vector  b  = 


A  A 

( kf)  *  . '  ■  *  bf  ) 
L0 


■l  II  b  y(n  -  i) |2. 


0  N=tg+1  *=0 


(2.7) 


Let  p^expOwj),  j  =  1,  ...,t0>  be  the  tQ  solutions  of 


A  °-  i 

B(z)  =  I  b.zJ  =  0, 
j=0  J 


(2.8) 


where  p.  >  0,  w.  e  [0,2ir),  j  =  1,  ...,tQ.  Then  we  use  w.'s  as  estimates  of 
J  J  ^  3 

I  c 

u>  •  S  . 

J 

Note  that  if  8^  satisfies 


st  =  BTt  l  I  I  b[t)y(n-£)|2,  t  =  0,l,...,T, 
z  n  rn=t+l  1=0  1 


(2.9) 


then  S.  is  the  smallest  eigenvalue  of  the  matrix 


f(t)  =  (^>),  t.™=  0,  1 . t, 


(2.10) 


'(t) 

and  bv  is  the  corresponding  unit  eigenvector,  where 
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1 

N-t 


N 


l  y(n  -  «.)y(n  -  m) , 
n=t+l 


i,  m  =  0, 1 ,  . 


t. 


(2.11) 


The  above  method  combines  the  procedure  of  parameter  estimation  with 
the  detection  procedure  of  the  parameter  number.  As  shown  in  Section  3, 

(tg*  0  1  tQ>)  is  a  strong  consistent  estimate  of  (tQ,  {w.,  j  <  t  ) 

2  ^ 

under  the  condition  (2.2).  Besides,  as  an  estimate  of  a  ,  S?  is  also 

l0 

strongly  consistent.  In  Section  4,  we  obtain  the  limiting  distributions 

of  S;  and  {to. ,  j  <  tn). 

Lq  J  u 
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3.  STRONG  CONSISTENCY  OF  THE  DETECTION  AND  ESTIMATION  PROCEDURES 
In  this  section,  we  establish  the  following. 


THEOREM  3.1.  Suppose  that  {w(n)}  is  an  i.i.d.  sequence  such  that  (2.2) 
holds.  Then  with  probability  one  for  N  large,  the  following  results  hold: 


(ii)  there  exist  a  unique  (tQ  +  l)xl  unit  vector  b  (up  to  a  complex 
factor  with  modular  one)  which  satisfies  (2.7),  and 

(iii)  for  appropriate  ordering, 

-  .  2 

oj  .  ->■  u  . ,  j  =  1  ,  . .  . ,  tn ,  S;  a  ,  as  N  ->■  ®. 

J  J  u  *-g 

In  other  words,  (tQ,  1^,  J  £  tQ},  S£  )  is  a  strongly  consistent  estimate  of 

2  ® 

(tg,  (Wj»  J  £  Tq  } ,  a  )  • 

To  prove  this  theorem,  we  need  the  following  lemma: 


LEMMA  3.1  (Petrov).  Let  {X  ,  n  >  1 }  be  a  sequence  of  independent  real 

n  ~  n  2  n 

random  variables  with  zero  means.  Write  s„  =  1  EX-  and  S  =  1  X..  If 

n  j=l  J  j=l  J 


lim  inf  s„/n  >  0 
n—  " 


and 

E | X.  12+6  <  K  <  •,  j  >  1 

J 

for  some  constants  K  and  i  >  o,  then 

lim  sup  S  / ( 2s^  log(s^))1//2  =  1,  a.s. 
n-x»  n  n  ri 


For  a  proof,  the  reader  is  referred  to  Petrov  (1975)  and  Stout  (1974,  p.274). 


Proof  of  Theorem  3.1. 


Under  the  model  (2.1),  we  have 


For  co j  f  9 
1 

N-v 


Thus 


By  condition 


i  N  _ 

-r  l  y(n-«,)y(n-m) 

l* _ 4-  i  T 


n=t+l 


N  l0 


nit  l  (  I  a .exp(i (n-a)w.)  +  w(n-t))(  l  a  -exp(-1  (n-m)io . )  +  w(n-m)) 
N  z  n=t+l  j=l  J  J  j=l  J  J 


t0  2  t0  ,  N 

l  I  a  .  |  exp(i(m-£)w.)  +  £  a  .a.  exp(i  (mu.- . )  )rrLF  >  exp(in(c 

j=l  J  J  j,k=l,j7kJIC  k  J  N_t  n=t+l 


-0  ,  N  _ _ _ 

+  [  a  .exp(i  (m-t)w.)  rr-i  J  exp(i  (n-m)u.)  w(n-m) 
j=l  J  J  T'n=t+1  J 


;o_  ,  N 

+  l  ajexp(i(m-0«j)  A-  l  exp(-1'n-A)u.)w(n-£) 

J=1  n=t+l  J 


+  KTT  1  w(n-s.)w(n-m) 


N 


N-t 


n=t+l 


^1  +  ^2N  +  ^3N  +  ^4N  +  ^5N  ( »  t,m=0,  l,...,t.  (3.2) 


j  f  k. 


l  exp(in(wj-oj|<) )  = 


t+1 


exp( i  ( t+1 )(-... -to,  ))  -  exp(i  (N+l ) (oo.-uj.  ))  , 

- J  K  _ J — L.  =  n(I) 


(N-t)  ( 1  -  exp(i(Wj.-W|<))) 


j  ^  k  ,  J  ,  k“l  ,  .  .  .  ,  tg  . 


J2N  = 


(3.3) 


(2.2)  and  Lemma  3.1, 
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By  the  law  of  the  iterated  logarithm  of  M-dependence  sequence, 


a2  +  0{. 


,  for  z  f  m, 


i ,  for  z  =  m. 


(3.5) 


flt)  *  3  Aim  +  I  |aj|2exp(i(m-t)»j). 

J  * 


z,  m  =  0,1,  . . . ,  t, 


(3.6) 


where  S„  is  the  Kronecker  sign. 
£m  3 


Using  (3.2)-(3.6), 


-  y(^  +  0(t 
m  m  4 


,  a.s. ,  z,  m  =  0,  1 ,  . . . ,  t. 


(3.7) 


Let  xjt}  >_  ...  >_  x[l]  be  the  eigenvalues  of  r^,  and  xj*^  >_ 
the  eigenvalues  of  r^.  Then 


>  x(t)  be 
•  1  At+1  De 


l  >  trr^r^L 

j-l  J  J  “ 


(see  von  Neumann  (1937)).  Hence, 


;(t)  _  x(t)  tr(f(t)  .  r^))2. 


I  (a  l'-A  <  tr(r< 
j=l  J  J 


(3.8) 


exp(-iw  )  ...  exp(-ios  )  ..  r 

n  1  »  A  =  di agLa, ,  a2,  . 

(t+l)xtQ  ...  . 


,  at  ], 
r0 


.  i 

exp(-itw-|)  ...  exp(-it«t  Jj 


(3.9) 
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Since  a^.  f  0,  j  =  1,  . ...  tg,  and  ^  t  wk  if  j  f  k,  it  is  easily  seen  that 


rank(QA)  =  min(t+l ,  tg) , 


(3.10) 


=  a2It+1  +  nAA*a*, 


(3.11) 


and  n*  denotes  the  transpose  of  the  complex  conjugate  of  fi.  Thus 


SM  >  o2 


At+1  >  °  f0r  1  "  V 


Xt+1  =  °2  for  t  -  V 


(3.12) 


_  ?(t) 


Hence,  by  (3.7)  and  (3.8),  noticing  that  St  =  X^j,  we  have 


ftl  2 

lim  St  =  x[+>  >  a  ,  a.s.  for  t  <  tQ, 


(3.13) 


|St-a 


,  a.s.  for  t  >  tr 


(3.14) 


Assume  that  t  <  tg.  Then  by  (3.13),  (3.14),  (2.6)  and  (2.5), 


lim(Rt  -  R. )  =  a2  -  <  0,  a.s. 

N-voo  XQ  X  11 


(3.15) 


Hence,  with  probability  one  for  N  large 


R*  <  R,.,  for  t  <  tr 
l0  z  1 


(3.16) 


Now  we  assume  that  t  >  tg.  Then  by  (3.14),  (2.6)  and  (2.5),  with 


probability  one  for  N  large, 


11 


RtQ  Rt  =  StQ  ■  St  ‘  ^  "  VCN 


.0(/STS!)  -  (t-t0)cN  <  o. 


(3.17) 


Hence,  with  probability  one  for  N  large. 


*N  =  V 


(3.18) 


which  establishes  (i) . 

To  prove  (ii)  and  (iii),  without  loss  of  generality,  we  can  assume 


tQ =  tQ.  By  (3.10)  and  (3.11)  with  t  =  tQ, 


<tn>  (tft) 


2  _  '"O'  '"O'  ,(t0} 

V1  to  ~  “  1 


(3.19) 


Thus,  the  equation 


ttn)  9 

(r  u  -  a- I  +1)b  =  0, 
V 


=  l 


(3.20) 


has  a  unique  root  b  =  (bg,b^,  . . . ,  b^.  )’  (up  to  a  complex  factor  with  modular 


one).  By  (2.7),  b  i‘s  a  unit  eigenvector  corresponding  to  the  smallest  eigen- 

a  +1  of  r  .  By  (3.7),  (3.8)  and  (3.19),  with  probability  one  for  N  large, 
L0 


j(V 

l0  1  z0  1 


(3.21) 


which  implies  that  the  equation 


.(tn)  „(tn) 

(r  +1)b  -  o. 


L0  '  "0 


(3.22) 


A  A 


has  a  unique  root  b  =  (bg,b-|,  ...,bt  )'  (up  to  a  complex  factor  with  modular 
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one),  and,  with  appropriate  choice  of  this  factor,  we  have 


lim  b  =  b,  a.s. 
N-~>  ~  ~ 


(3.23) 


Now,  choose  b  *  (bn,  b. ,  ...,  b.  )'  such  that  ||b||  =  1  and 
u  1  **0 

A  t()  t 

B(z)  =  l  b.z*  =  0  (3.24) 

£=0 

has  tg  roots  exp(iu^) ,  expO'o^) ,  ...,  exp(iujt  ).  Then  b  is  the  unique  root 
of  equation  (3.20)  (up  to  a  complex  factor  with  modular  one)  and  vice-versa. 
Using  (3.23), 


lim  B(z)  =  B(z),  a.s.  (3.25) 

N-»°° 


Now,  use  the  definition  of  p.exp(id).),  j  =  1,  ...,tQ,  for  appropriate 
ordering, 

lim  p.exp(iw.)  =  exp(iw.),  a.s.,  for  j  =  1,  ...,tn, 

N-h»  j  j  j 

which  implies 


1  tm  at.  -  w . ,  a.s.,  for  j  l,...,tp.« 

N+oo  J  J  u 


(3.26) 


Using  (3.14),  we  establish  parts  (ii)  and  (iii)  of  the  theorem. 


Remark.  The  EVLP  method  can  be  easily  generalized  to  EVFBLP  along  a 
similar  line.  The  anoloques  of  the  results  qiven  in  this  section  and  the 
next  section  can  be  proved  step  by  step  as  the  proof  given  in  these  two 
sections.  Also,  we  can  expect  EVFBLP  to  give  more  accuracy  when  SNR  is 
high  and  the  sample  size  is  small. 


<:■  ,'v'vVv 
-■^.-1 


<■  H.l 


'i  t«j  ,»<  tu  «<*  ■nvuvt'.iiMt'.u'.M  «♦» 


.  jrnifmnurmrm.'i 
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4.  LIMITING  DISTRIBUTION  OF  S;  AND  (£.,  j  <  tn) 

r0  J 

Since  tg  -*■  tg,  a.s.,  as  N  -*■  »,  we  can  use  tg  instead  of  tg  when  we 
consider  the  limiting  properties  of  various  estimates  involving  t  For 
further  simplicity  of  notation,  we  will  omit  the  subscript  0  of  tg  and  sim¬ 
ply  use  t  for  tg.  Also,  we  will  keep  all  other  notation  defined  in  pre¬ 
vious  sections.  In  particular,  the  matrix  a  in  (3.9)  is  a  (t  +  l)xt  matrix. 

Throughout  this  section,  we  assume  that  (w(n)}  is  a  sequence  of  i.i.d. 
complex  random  variables  such  that 

Ew(l )  =  0,  E(Rw(l))2  -  E(Iw(l))2  =  \o2, 

E ( Rw ( 1 ) Iw( 1 ) )  =  0,  and  Var(|w(l)|2)  =  ao^  with  a  >  0.  (4.1) 

Here,  Rw(l)  and  Iw(l)  denote  the  real  and  imaginary  parts  of  w(l)  respec¬ 
tively. 

LEMMA  4.1.  Suppose  that  condition  (4.1)  is  satisfied.  Then 

1  N  ■n 

-  y  exp(-i(n  -  Ow.Mn  -  *,)  — ►v.,  j  =  l,2,...,t,  £  =  0,1,  ...,t4 

/fPt  n=t+l  J  J 

N 

--  ^  £  (|w(n-£)|2-o2)  ■ — *■  un ,  £  =  0,  1 ,  . . . ,  t, 

/N-t  n=t+l  J 

i  N  _  v 

-  I  w(n  -  £)w(n  -  m)  — *•  u  ,  if  0  <  m  <  £  <_  t. 

n=t+l  x'"m 

Here  v.'s  and  u.'s  are  independent  of  each  other  and 
J  J 

(i)  Vj  -  Nc(0,  o2),  j  =  1,  ...,  t, 

(ii)  Ug  -  Nr(0,  ao4), 

( i i i )  u  .  -  N  (0,  o4) ,  j  =  1 ,  . . , ,  t .  (4.2) 
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We  have  the  following  result: 

!  THEOREM  4.1.  Suppose  that  (w(n)}  is  a  i.i.d.  sequence  of  complex 

j 

!  random  variables  such  that  (4.1)  holds.  Then  we  have 

i 

j  ?  =  b*Ub,  (4.5) 

>  and 

» 

!  Tn  +  iAN  G’'I(AA*)“1(ft*fl)"1Q*Ub.  (4.6) 

i 

Proof.  Put  V  =  diag[v-j ,  . . . ,  v^.]  and 

H  =  /Nht(f(t)  -  r(t))  =  (hm).  (4.7) 

By  (3.2),  (3.3)  and  Lemma  4.1,  we  have 


Since  St 


H  -2+  H  *  (h£m)  =  q(AV  +  AV)n*  +  U. 

~  ( t) 

the  smallest  eigenvalue  of  r'  ,  we  have 
0  =  det[r(t)  -StIt+1] 

=  det[r(t)  -  o2It+1  +  — - —  H  -  (S.  -  o2)It+1] 

/fPt  z  1  1 

=  det[ftAA*n*  +  — —  H  -  (S.  -  a2)It.,]. 

/Pt  z  z 


(4.8) 


(4.9) 


Let  Q  be  a  unitary  matrix  such  that 

Q*fiAA*fi*Q  =  d i a g [ C i ,  . . . ,  0],  C1  >  •••  >  Ct  >  0. 

Note  that  the  last  column  of  Q  is  the  eigenvector  of  corresponding  to 
2  (t) 

the  eigenvalue  a  of  rv  .  We  can  choose  this  column  as  b.  Now  we  have 


I  ISJ'i  I*,  it,  1 


0  =  det 


5r(st"°  ) 


V*st"°  )  i 

z  z  +  — ^  Q*HQ 

/c  2,  ATt 

'< st"°  ]  J 
/ 


(4.10) 


Since  H  H,  by  Skorohod's  representation  theorem  (see  [7]), 
we  can  assume  H  H,  a.s.  as  N  -*•  «.  Multiply  by  (N-t)^  the  last 
row  and  the  last  column  of  the  matrix  in  {  }  of  (4.10)  respectively.  By 


(4.10),  $t  -*•  o  ,  a.s.  and  noting  that  fi*b  =  0,  we  have 


-*•  c,  =  b*Hb  =  b*Ub,  a.s. 


(4.11) 


Note  that  (4.11)  reveals  only  that  there  exist  some  versions  of  c, 

H  and  U  which  have  the  same  distributions  as  c ,  H  and  U  respectively 
such  that  (4.11)  holds.  From  this  we  only  get  (4.5).  The  principle  of 
this  statement  also  applies  to  the  following  proof  of  (4.6)  and  so  on. 

Since  (r^  -  StIt+1  )b  =  0  and  (r^  -  ct2I t+1  )b  =  0  with  the  choice  of 
b  such  that  b  -*•  b,  a.s.,  we  have 


^t(r(t)  "Vt+rl? 

^Pt(r(t).02r  +  _Lh — L_?  i  )(b+-i—n  ) 

t+1  /FPt  ATt  N  t+1  '  /fPt~N 


(r^-02 


a  Wjn  +  1 1+ 1 +  (H-cNit+1)(b-b), 


where 


^b-b)  *  (nN0,...,nNt)’. 


(4.12) 


ite  nN  =  such  that  ni1  ^  €  u(n)  and  ]_  u(n) ,  where 


u(ft)  denote  the  space  spanned  by  all  column  vectors  of  n.  Then 

for  some  txl  complex  vector  s^.  Since  H  -*■  H,  a.s.,  -*■  t,  -  b*Ub,  a.s. 

and  b  b,  a.s. ,  we  have 

(H  -  C^t+l )(’b  -  b)  -►  0,  a.s.  as  N  -*■  «.  (4.13) 

Note  that 

(r^  -o2It+1)riN  =  nAA*n*jN  =  aAA*n*n^  =  nAA*n*ngN  (4.14) 

and 

(H  -  ^N!t+1 "  (H*cIt+l)-  =  U-  " 

=  (It+1  -bb*)Ub,  a.s.  (4.15) 

By  (4.12)-(4.15) ,  and  It+1  -  bb*  =  n(n*n)_1n*,  we  get 

«N ->  -(n^r1(AA*r1(a*n)”',n*Ub,  a.s. 
which  implies  that 

-*■  -fi(Q*n)_1(AA*)_1  (n*n)-1n*Ub,  a.s.  (4.16) 

Hence 

a*nN  =  *  -(AA*r1(fi*n)'Vub,  a.s.  (4.17) 

Finally,  let  us  consider  the  limiting  distribution  of  A^.  Since 

A 

B(ojexp(iij))  =.  B(exp(iu.))  =  0,  we  have  for  j  =  l,...,t. 


i  • 
I 

i 


| 
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0=^1  b  p!'exp(Uw.) 
£=0  *  J  J 

t 


where  p.*  e  [l,p.]  or  [p.,1]  and  w  •*  e  [w.,w.]  or  [w.,w.].  We  have 
J  J  J  J  J  J  J  J 

proved  that  b  ■+  b,  a.s.,  p .  -*•  1 ,  a.s.  and  u>.  to.,  a.s.  for  appropriate 

\J  s)  \J 

ordering  and  for  j  =  1,  t,  so  that,  with  probability  one,  we  have 

fi*nN  +  G<In  +  1V  +  0  ^  1 1  1 1 +  II  II  )  =  9  (4.18) 

for  large  N.  Since  exp(ioj.)  is  a  simple  root  of  the  multinomial  B(z),  we 

J 

have  D(exp(iw.))  ?  0  for  j  =  1,2,  ...,t,  and  6  is  nonsingular.  From  this, 

J 

and  (4 . 1 7)- (4 . 1 8) ,  it  follows  that 

IN  +  iAN  G"1(AA*)”1(n*n)“1a*Ub,  a.s.  (4.19) 


As  mentioned  above,  this  only  gives  the  corresponding  weak  convergence 
described  in  (4.6).  This  completes  the  proof  of  Theorem  4.1. 


i 

» 

I 

» 


5.  FURTHER  DISCUSSION  ON  MODIFIED  FBLP  METHOD 


Consider  the  model  (2.1)  and  suppose  that  (4.1)  holds.  In  this  section 

we  discuss  the  problem  of  estimating  the  frequencies  when  tg  is  known.  It 

is  emphasized  that  some  of  the  notations  of  this  section  may  be  different 

from  those  of  the  previous  sections. 

Nutall  (1976)  and  Ulrych  and  Clayton  (1976)  developed  the  method  of 

forward-backward  linear  prediction  (FBLP),  which  works  well  when  the  signal  - 

to-noise  ratio  (SNR)  is  sufficiently  high.  However,  it  becomes  false  in 

the  case  of  low  SNR.  Tufts  and  Kumaresan  (1982)  proposed  a  modified  FBLP 

(MFBLP)  method.  They  showed  by  simulation,  that  the  estimation  of  the 

frequencies  has  been  greatly  improved  by  MFBLP  when  the  SNR  is  relatively 

low.  In  this  section  we  shall  give  a  theoretical  analysis  to  show  why  the 

MFBLP  improves  when  SNR  is  low  and  the  sample  size  is  large.  Also,  we 

shall  suggest  a  further  modification  to  this  method. 

First,  let  us  describe  the  FBLP  procedure.  The  reader  should  note 

that  it  is  assumed  that  the  true  number  of  signals  is  known  in  this  procedure. 

-ft) 

The  linear  prediction  coefficient  vector  bv  (t>_tQ)  is  defined  as  the 
one  which  minimizes 

N  t  N-t  t 

l  |y(n)  +  l  b^y(n-il)  !2  +  l  |y*(n)  +  f  bjt^y*(n+i)  |2,  (5.1) 

n=t+l  £=1  1  n=l  £=1  * 

where  b^  ,  i  =  1 ,  2,  . . . ,  t,  are  the  components  of  b'  . 

In  case  the  solution  is  not  unique,  we  select  the  one  which  also  mini¬ 
mizes 

I  |b‘l)|2.  (5-2) 

£=1  1 

A(t) 

In  light  of  this,  we  call  the  solution  bv  '  least  normed  prediction  (LNP) 


coefficient  vector. 


We  can  construct  the  transfer  function  of  the  prediction  error  filter 


B(t)  (z)  =  1+1  b;wz'\ 


(t)  -i 


(5.3) 


£=  1 


Let  pSt)exp(i^t)),  p\(t)  >  0,  ^t)  e  [0,2*),  £  =  l,2,...,t0  be  the  first 


t_  roots  of  6^(z)  which  are  closest  to  the  unit  circle  in  the  complex 


plane  among  the  t  roots,  and  take  £  =  1,  2,  ...  t«  as  the  estimates  of 


1  c 

W£  S* 


To  introduce  the  MFCLP  method,  write 


and 


Y  = 


y(t-D. 

y(t). 


•  •  •  » 

•  •  •  9 

•  •  •  9 


y(N-2), 
y*(3) , 
y*(4) , 


•  •  9 

•  •  9 
■  •  •  9 


y*(N-t+l), 


y*(N-t+2), 
;(t)  ^ 


;(t)  . 


/ 


'1 


h  =  (y(t+l ) ,  y(t+2) ,  y(N),  y*(U.  ■•••  y*(N-t)) 


It  is  well-known  that  the  LNP  is 


=  -(Y*Y)+Y*h 


(5.4) 


where  (Y*Y)  denotes  the  Moore-Penrose  inverse  of  (Y*Y). 

Now,  let  (2N-2t)”^Y*Y  have  the  following  spectral  decomposition 


(2N-2t)~  Y*Y  =  T  i.u.u* 
j=l  J'J"J 


(5.5) 


As  shown  in  the  sequel,  with  probability  one  for  large  N,  the  above  matrix 
is  positive  definite.  Hence,  the  FBLP  gives  the  LNP  as  follows: 


i  q1  <«);,. 

j=l  J  'J~  -J 


where 


=  -(2N-2t)_1Y*h. 


(5.6) 


(5.7) 


If  we  use 


0 

b*  =  l  V(uSy)u, 


J  -J~  ~J 


(5.8) 


instead  of  b^  in  the  above  FBLP  procedure,  we  get  the  MFBLP  method. 
In  this  section,  we  propose  to  use 


l  (X--A  )'  (u]y)U-, 

j=1  J  t  'J'  'J 


t  >  t0’ 


(5.9) 


instead  of  b*  in  the  MFBLP  procedure.  We  will  establish  the  strong  consis¬ 
tency  of  the  estimates  of  the  frequencies  for  this  method  under  the  condi¬ 
tion  (4.1). 

Now,  we  consider  the  case  without  noise.  To  distinguish  this  case  from 


the  case  with  noise,  we  drop  out  the  superscript  'V  in  all  notations 
Wri  te 


and 


Then 


and 


Q  = 


a1exp(itcj-|) , 

. . . ,  af  exp(itan  ) 

L0  r0 

a1exp(i(t+l)w1), 

•  •  •  > 

. . . ,  a .  exp(i  (t+1  )oj . 

z0  z0 

a.expfifN-l)^), 

....  a.  exp ( i ( N- 1 W 

z0 

a1exp(-2ioJl), 

...,  a.  exp(-2iwt  ) 

a'1exp(-3ito1 ) , 

•  •  •  9 

...,  a.  exp(-3iw,  ) 
r0  r0 

•  •  •  9  •••» 

a’1exp(-i  (N-t+Dcj^ ) , 

. . . ,  a.  exp(-i (N-t+1 ) 

Q  = 


expHwj) 

exp(-iw2) 


•  •  • 

•  •  • 


1  exp(-iw.  ) 


Y  =  Qft 


exp(-i  (t-1  Jw-j ) 
exp(-i(t-l)tD2) 


I 

PXp(-l'  (t-1  )  UK  ) 

0  / 

/ 


b  -  Q0g» 


where 


20  =  (exp(ioj1),  ....  exp(iwt^) ) ' . 


Substituting  (5.10),  (5.11)  into  (5.4),  we  get 

=  - ( q*Q*Qq) +&*Q*QnQ . 


-y.j 


LA. 


<< 


(2N  -  2t)  'Y*h  =  fi*An0  +  0(/^log  logN),  a.s. 


(5.24) 


Note  that  tr1'  is  independent  of  N  in  view  of  (5.13).  Hence,  by  (5.12)  and 
(5.22)  we  get 


=  -(ft*Afl)+n*Ag0 


=  (ft*Afi)+e. 


(5.25) 


Here  we  have 


§  =  -fi*AfiQ. 


(5.26) 


Note  that  the  rank  of  the  txt  matrix  is  tQ.  So  the  matrix  ft*Aft  has 


the  following  spectral  decomposition 


U 

n*Afi  =  l  C..-V.V  * 

j=l  J  J  J 


(5.27) 


where  >  ^  >  ...  ^  are  the  non-zero  eigenvalues  of  n*An,  and  v.'s  are 
1  c  ‘•n  '3 


the  corresponding  eigenvectors. 


Applying  Lemma  2.1  in  Bai-Krishnaiah-Zhao  (1985),  by  (5.23),  we  obtain 


aj  =  Cj  +  +  °(/^log  log  N) ,  a.s.,  for  j  <  tQ, 


Xj  =  o'  +  0(  /jq-log  log  N), 


a.s. ,  for  j  =  tQ+l ,  . . . ,  t.  (5.28) 


By  (5.5),  (5.23),  (5.27)  and  (5.28),  we  get 


v.  2  -  -  u 

y  (xH  -  o  )u  yt  =  y  c.v.v*  + 

j=l  J  J  J  j=l  J 


-) »  a.s. 


+  0(' 


Hog  log  N 
N 


).  a.s. 


(5.29) 


wjtowkw  jpi*v*wvjiiuuui  «uu 


26  . 


and 


+  “Vtlr^  +  0(/i£MraJ1  >•  a-s- 

«J_ln  1 


(5.30) 


Take 

¥j»  i  =  V1 

i  •  •  •  )  t  dS  the 

rest  unit  eigenvectors  of 

such  that 

v*v„ 

~j~* 

=  0  for  j  f 

t.  j,  £■ 

1  »  •  •  • 

t.  Put 

G1  " 

(ur  ... 

5  U.  ), 
** 

U2 

=  ^t«+l *  ** *  ’ 

5t>- 

txtQ 

U 

tx(t-tQ) 

0 

V1  ' 

(v-j.  ... 

»V.  ), 

^  tri 

h 

^-tn+l *  •*•* 

txt0 

0 

tx(t-t0) 

0 

and 

* 

U2  " 

V1  S 

+  V2 

G2N* 

tx(t-tQ) 

txt0 

tx(t 

"t0) 

where 

G1N  and  G2N 

are  tQ x 

(t-tQ) 

and  (t-tg) 

x  (t-tg)  matrices  respectively.  | 

By  (5.29), 


-) ,  a.s.  for  k  =  tg+1 ,  .. . ,  t. 


Thus 


Gfv^oi/SEal),  ,.s. 

By  (5.30)  and  (5.31), 


(5.31) 


■'Ml  ■  °VtQ  *  0(/^S).  a.s. 


(5.32) 


Thus 


Similar  to  (5.31),  we  have 


v*u2  = 


,  a.s. 


which  implies  that 


°1N  ■  °<, 


-),  a.s. 


From  the  above  argument,  it  follows  that 


l  u.iit  =  G?Ui  =  V?V*  + 
j=t0+rj~J  c  c  c 

*  .1  Vj-Yj  +  c 
J=tn+1  J  J 


i  °  U.Ut 


j=V] 


?  l  v.v*  +  0(c" 

j=tn+rJ'J 


)  >  d .  S . 


,  a.s. 


!  :-i 


y  V •  o  I  y-vt  +  ou' 

J-tn+i  J  V‘ 


,  a.s. 


(5.34) 


(5.35) 


(5.36) 


(5.37) 


(5.38) 


Here  we  write  the  factor  <j  in  the  remainder  in  order  to  compare  the  above 
methods  in  the  sequel.  Using  similar  argument,  we  can  prove  that 


V,~  -  *o 

V  ’it  n*  =  V 


l  u.ut  =  l  U-+a  )  v.vt  +  o( 
j  =  l  J  j=i  J  J  J 


to  .  .  .  *0  , 


,  a.s. 


(5.39) 


,  a.s. 


(5.40) 


» * . ' .  ■ 


By  (5.7),  (5.24)  and  (5.26), 


Y  =  3  +  0( 


,  a.s. 


(5.41) 


Note  that  by  (5.25)  and  (5.27),  we  see  that 


■*in\  +  o  = 


l  cV(v*e)v.  =  (n*An)  b  =  b 
j=!  J  ~J~  ~J  ~  ** 


(5.42) 


By  (5.27) ,  we  have 


V|ft*An  =0,  and  V|Q*A  =  0. 


From  this  and  B  =  -n*l\nn,  it  follows  that 
-  -U 


(5.43) 


1  (v^b)v  .  =  0. 
,=t  +1  ~3~  -3  - 


(5.44) 


Thus,  by  (5.6),  (5.8),  (5.9),  (5.38)-(5.42)  and  (5.44),  we  have 


b^  =  l  U.j+cj^)  ^(vts)v.  +  0(a  ^ 

•  ^1  J  ~J  ~  -J  » 

J 


b*  =  l  U-  +  a2)_1(vtB)v.  +  0( 
j  =  !  J  ~J~ 


b  =  b^  +  0( 


~)  »  d  .  S  . 


L) ,  a.s. 


(5.45) 


(5.46) 


,  a.s. 


(5.47) 


Put  b  =  (b, ,  . . . ,  b. ) 1 ,  B(z)  =  1  +  l  b  z‘J  and  let  ojexp(i^.) , 

1  z  z=l  1  3  3 

A  /S  /V 

j  =  1  f  2,  . . .  f  t,  be  the  roots  of  B(z) ,  where  i  P2  1  >  0*  As 

pointed  out  earlier,  w  ^ ,  j  =  1,  ...,tg,  are  the  roots  of  B^(z).  Kumaresan 

(1982)  has  shown  that  the  particular  choice  of  b^  places  the  t  -  tg  extraneous 


# 

‘,1  ,\> 


WwlvIW 
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zeros  of  the  prediction-error  filter  transfer  function  B^(z)  inside 
of  the  unit  circle.  Thus,  by  (5.47)  we  obtain  that  for  appropriate  order- 

A  A 

1  HQ  Of  W’ s  •  •  •  »  (i)i,  * 

1  ZQ 

“j  8  “j  +  0(/^f^)>  a.s.  for  j  =  1 ,  . . . ,  tQ,  (5.48) 

which  implies  the  strong  consistency  of  the  third  method  by  us. 

By  (5.45)  and  (5.46),  we  see  that  the  FBLP  procedure  and  the  MFBLP 
procedure  have  the  same  asymptotic  behavior. 

Denote  by  B(z)  the  transfer  function  constructed  according  to  (5.3) 

■'(tl  0  2-1 

with  b'  '  replaced  by  £  (c .  + a  )  (v*b)v  . .  In  the  same  way,  we  can  con- 

j=l  J  *J~ 

*  *  A  /  f  \ 

struct  B*(z)  using  b*  instead  of  bv  '  in  (5.3).  By  (5.46)  we  know  that 

B*(exp(iUj))  =  B(exp(iuj))  +  0(/^ffi-N),  a.s.  (5.49) 

for  j  =  1 ,  . . . ,  tQ.  In  general ,  B(exp(iw.) )  +  B^  (exp(iw . ) )  =  0  for 
j*l,  ...,tg  in  view  of  (5.42),  so  that  MFBLP  method  is  not  consistent. 

But  we  know  that  B(z)  reduces  to  Bv  ' (z)  when  a  =  0.  Thus  when  o  is  small 
enough  (with  respect  to  >.,  j  =  1 ,  . . . ,  tQ) ,  B(exp(iw.))  =  0  and  B*(exp(iw  ■) ) 

2  0  for  j  =  1 ,  2,  . . . ,  tg. 

In  the  same  way, 

B^(exp(iuj.))  =  B(exp(iw,))  +  0(o~2  ) ,  a.s.  (5.50) 

\J 

2 

for  j  =  1 ,  . . . ,  tg.  When  a  is  small ,  the  main  term  of  the  above  expression 

2 

B(exp(ioj.))  -  0  for  j  <_  tg.  This  means  that  when  a  is  small  and  fixed, 

3  — 

the  FBLP  procedure  can  estimate  the  true  frequencies  well  in  the  large 
sample  case.  However,  when  N  is  not  large,  the  remainder  in  (5.50)  would 


3i 


bring  much  random  effect  on  the  zeros  of  the  prediction-error  filter 

transfer  function  B^(z).  As  shown  by  the  simulation  results  given  in 

Fig.  4  and  Fig.  5  in  Tufts  and  Kumaresan's  paper  (1982),  the  main  random 

effect  is  imposed  on  the  t  -  tn  extraneous  zeros  of  B^(z).  If  we  drop 
t  u 

l  A.  (u?y)u.  when  we  construct  the  transfer  function,  as  done  by  Tufts 
j=tQ+l  J  ~J~  ~J 

and  Kumaresan,  the  estimates  of  the  true  frequencies  certainly  can  be  im¬ 
proved.  This  is  why  the  MFSLP  procedure  is  better  than  the  FBLP  procedure, 
as  shown  by  the  above  simulation  results. 

Thus,  we  establish  the  following. 

THEOREM  5.1.  Suppose  that  (4.1)  holds  in  the  model  (2.1).  For  the 
third  procedure  proposed  by  us  in  this  section,  and  for  appropriate  order¬ 
ing  of  w, ,  ....  ui.  ,  we  have 
i  tQ 


"j  ■  “j  *  Ol/SSI),  a.s. 


for  j  =  1 ,  . . . ,  t. 
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